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Direct  numerical  simulation  (DNS)  results  of  turbulent  MILD  premixed  and  conventional  (undiluted) 
premixed  combustion  have  been  investigated  to  shed  light  on  the  physical  aspects  of  reaction  zones 
and  their  morphology  in  MILD  combustion.  Results  of  a  premixed  case  are  used  for  comparative  analyses. 
The  analyses  show  that  the  regions  with  strong  chemical  activity  in  MILD  combustion  are  distributed 
over  a  substantial  portion  of  the  computational  domain  unlike  in  the  premixed  case  where  these  regions 
are  confined  to  a  small  portion  of  the  domain.  Also,  interactions  of  reaction  zones  are  observed  in  MILD 
combustion  with  their  spatial  extent  increasing  with  dilution  level.  These  interactions  give  an  appearance 
of  distributed  combustion  for  MILD  conditions.  The  morphology  of  these  reaction  zones  is  investigated 
using  the  Minkowski  functionals  and  shapefinders  commonly  employed  in  cosmology.  Predominant 
sheet-like  structures  are  observed  for  the  premixed  combustion  case  whereas  a  pancake-like  structure 
is  observed  as  the  most  probable  shape  for  the  MILD  cases.  Spatial  and  statistical  analyses  of  various 
fluxes  involved  in  a  progress  variable  transport  equation  are  conducted  to  study  autoignitive  or 
propagative  characteristics  of  MILD  reaction  zones.  The  results  suggest  that  there  are  local  regions  with 
autoignition,  propagating-flames,  and  their  coexistence  for  the  conditions  considered  in  this  study. 
Typically,  reaction  dominated  or  ignition  front  and  propagating-flame  dominated  regions  are  entangled 
for  high  dilution  cases.  Scalar  gradient  plays  a  strong  role  on  whether  reaction  or  propagating-flame 
dominated  activities  are  favoured  locally. 
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1.  Introduction 

Moderate  and  Intense  Low-oxygen  Dilution  (MILD)  combustion 
has  the  potential  to  improve  combustion  efficiency  and  reduce 
pollutant  emissions  simultaneously  [1-3  .  This  type  of  combustion 
is  associated  with  intense  preheating  and  dilution,  where  (i)  the 
elevated  reactant  temperature,  Tr,  is  higher  than  the  autoignition 
temperature,  Tign,  of  a  given  fuel  mixture  and  (ii)  the  temperature 
increase  during  combustion,  AT  =  TP  -  Tr ,  is  lower  than  Tign.  The 
product  temperature,  Tp,  is  very  low  compared  to  that  in  conven¬ 
tional  combustion  even  with  the  elevated  reactant  temperature 
because  the  reactant  mixture  is  diluted  with  a  large  amount  of 
exhaust  gas  so  that  the  oxygen  level  in  the  reactant  mixture  is  typ¬ 
ically  low,  from  2  to  5%  by  volume  [3  .  MILD  combustion  is  clearly 
distinguished  from  conventional  premixed  combustion  since  Tr 
and  AT  are  quite  different  in  the  two  modes  of  combustion  [3  . 
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Combustion  under  MILD  conditions  has  many  advantages  [1-3] 
as  summarised  briefly  below.  Firstly,  combustion  efficiency  is 
enhanced  by  preheating  the  reactants  using  the  heat  recovered  from 
the  exhaust.  Secondly,  NO  formation  is  suppressed  significantly 
because  of  intense  dilution  of  reactants  resulting  in  low  oxygen  level 
and  flame  temperature.  It  typically  takes  a  few  seconds  to  produce 
substantial  levels  of  thermal  NO  at  around  1900  I<  but  this  time 
reduces  to  a  few  milliseconds  at  2300  K  [1,2  .  The  maximum 
temperature  in  MILD  combustion  is  typically  less  than  1900  K 
[1,3  .  Thirdly,  combustion  conditions  having  large  Tr  and  small  AT 
help  to  suppress  combustion  noise  and  instabilities  [1,2  .  MILD 
combustion  has  been  achieved  in  various  configurations  in  previous 
studies  [1-10]. 

Previous  studies  employing  direct  photographs  [6,7,9]  and 
laser  thermometry  [4,5,7,11]  suggest  a  uniform  and  distributed 
combustion  under  MILD  conditions,  often  referred  to  as  flameless 
combustion.  On  the  contrary,  the  Probability  Density  Function 
(PDF)  of  temperature  reported  in  [8  suggests  the  existence  of  thin 
reaction  zones.  Planar  Laser  Induced  Fluorescence  (PLIF)  images  of 
OH  radicals  also  show  fairly  thin  reaction  zones  [4,5,7,11].  A 
recent  study  using  direct  simulation  data  of  MILD  combustion  of 
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a  methane-air  mixture  shows  the  existence  of  thin  reaction  zones 
which  interact  frequently  leading  to  the  appearance  of  distributed 
combustion  [12  .  Thus,  it  is  imperative  to  analyse  reaction  zone 
characteristics  to  understand  whether  MILD  reaction  zones  are 
thin  or  not  since  this  understanding  is  important  to  develop 
appropriate  modelling  for  turbulent  MILD  combustion. 

Autoignition  is  also  expected  to  play  a  role  in  MILD  combustion 
because  Tr  >  Tign  as  noted  in  previous  studies  [1-3].  The  role  of 
autoignition  on  flame  stabilisation  has  been  investigated  using 
the  Jet  in  Hot  Coflow  (JHC)  burner  [5,13-15]  and  in  other  configu¬ 
rations  involving  diluted  and/or  hot  reactants  and  coflow  16,17]. 
The  PLIF  images  of  CH20  in  5]  suggest  that  the  hot  coflow  initiates 
ignition  of  the  MILD  mixture  after  localised  extinction  is  caused  by 
the  entrained  cold  surrounding  air  in  the  JHC  configuration. 
However,  PLIF  images  of  OH  [4,5,7,11]  and  temperature  PDFs  [8] 
suggest  that  the  MILD  combustion  in  the  JHC  burner  involves 
propagating  flames.  Thus,  an  interplay  between  autoignition  and 
propagating  flame  can  be  expected  in  MILD  combustion  [3,18]. 

These  contradicting  views  on  MILD  combustion  involving  either 
thin  reaction  zones  or  distributed  combustion,  and  either  autoigni¬ 
tion  or  propagating  flames,  are  the  primary  questions  to  be  tackled 
in  the  present  study.  Direct  numerical  simulation  (DNS)  data  are 
ideally  suited  for  this  purpose  and  a  detailed  analysis  of  MILD  com¬ 
bustion  DNS  data  is  conducted  to  answer  these  questions.  The 
morphology  of  turbulent  MILD  reaction  zones  is  studied  using 
the  Minkowski  functionals  to  be  described  in  Section  3.2.  The 
autoignitive  or  propagative  characteristics  of  these  reaction  zones 
are  analysed  using  the  balance  of  various  fluxes  in  the  transport 
equation  of  a  reaction  progress  variable  based  on  temperature. 

This  paper  is  organised  as  follows.  The  DNS  methodology,  con¬ 
struction  of  initial  fields,  and  combustion  conditions  are  briefly 
described  in  Section  2.  The  results  are  presented  in  Section  3.  The 
PDFs  of  reaction  progress  variable  in  MILD  combustion  are 
presented  in  Section  3.1.  The  analysis  using  Minkowski  functionals 
is  discussed  in  Section  3.2.1  and  the  flux-balance  analysis  is 
presented  in  Section  3.3.  The  conclusions  are  summarised  in  the 
final  section. 

2.  DNS  of  MILD  combustion  and  turbulent  premixed  flames 

Combustion  of  a  lean  methane-air  mixture  will  be  considered 
under  both  MILD  and  conventional  premixed  conditions.  In  this 
paper,  MILD  premixed  combustion  will  be  referred  to  as  “MILD 
combustion”  while  conventional  undiluted  premixed  combustion 
will  be  referred  to  as  “premixed  combustion”.  The  conditions  for 
the  MILD  combustion  are  created  using  the  concept  of  exhaust 
gas  recirculation  (EGR)  and  thus  the  reactant  mixture  is  made  of 
unburnt  and  burnt  gases  distributed  randomly  in  space  and  time. 
The  construction  of  this  field  is  described  in  Section  2.2.  Direct 
simulation  of  a  complete  EGR  system  is  beyond  the  reach  of 
current  computational  capability  and  thus  a  two-stage  strategy  is 
followed  in  this  study.  Before  describing  this  strategy  in  detail  in 
Section  2.2,  the  numerical  methods  are  discussed  first. 

2.2.  Numerical  method  and  configuration 

The  numerical  code  SENGA2  [19  employed  for  this  study  solves 
fully  compressible  transport  equations  for  mass,  momentum,  total 
energy  and  mass  fractions  of  various  scalars  involved  in  the 
combustion  chemistry.  The  transport  properties  depend  on  tem¬ 
perature  and  a  detailed  chemical  kinetics  scheme  can  be  used  for 
combustion  chemistry.  This  code  has  been  used  in  several  earlier 
studies  on  turbulent  premixed  flames  [20,21]  involving  one-step 
or  complex  chemistry.  Here,  a  skeletal  mechanism  involving  16 
species  and  25  elementary  reactions  [22]  is  used.  Binary  Fickian 


diffusion  is  used  with  constant  Lewis  numbers  for  each  species. 
The  reference  autoignition  delay  times  computed  using  this 
diffusion  model  are  comparable  to  those  obtained  using  a 
multi-component  diffusion  model  for  MILD  conditions  [15  .  Also, 
the  use  of  multi-step  chemistry  and  non-unity  Lewis  numbers  is 
important  for  MILD  combustion  because  of  the  role  played  by  rad¬ 
icals  and  intermediate  species  present  in  the  diluent.  The  values  of 
burnt  side  temperature,  laminar  flame  speed,  thermal  thickness, 
and  reference  autoignition  delay  times  computed  using  the  skele¬ 
tal  mechanism  compared  well  with  those  obtained  using  GRI3.0 
mechanism  for  the  MILD  conditions  of  this  study. 

The  spatial  derivatives  are  approximated  on  a  uniform 
numerical  grid  using  a  tenth  order  central  difference  scheme  which 
gradually  reduces  to  a  fourth  order  scheme  near  boundaries.  The 
time  integration  is  achieved  using  a  third  order  Runge-Kutta 
scheme.  The  numerical  stability  of  these  schemes  is  maintained 
by  using  sufficiently  small  time  steps  (At  ^  1  x  1(T8  s),  which  are 
dictated  by  the  acoustic  CFL  condition.  The  computational  domain 
is  cubic,  and  inflow  and  non-reflecting  outflow  boundaries  in  the  x 
direction  are  specified  using  the  Navier-Stokes  Characteristic 
Boundary  Condition  23]  while  periodic  conditions  are  specified 
in  the  y  and  z-directions.  A  similar  computational  domain  with 
these  boundary  conditions  is  used  for  a  conventional  turbulent 
lean  premixed  flame,  which  is  used  for  comparative  analyses. 

The  mixture  fields  for  MILD  combustion  contain  pockets  of 
exhaust  and  fresh  gases.  Such  mixture  fields  are  generated  care¬ 
fully  in  preprocessing  steps  to  be  discussed  in  Section  2.2.  These 
preprocessed  fields  of  scalar  mass  fractions,  T/[x(t),y,z],  velocity, 
u[x(t),y,z],  and  temperature,  T\x(t),y,z]  are  fed  at  an  average  veloc¬ 
ity  of  Uin  as  inflowing  fields,  Yfix  =  0,y,z,  t),u(x  =  0,y,z,  t)  and 
T(x  =  0,y,z,t)  through  the  inflow  boundary  located  at  x  =  0  of 
the  computational  domain.  The  symbol  x(t)  denotes  the  x  location 
of  a  scanning  plane  at  time  t  moving  at  a  velocity  of  Uin  through  the 
preprocessed  fields.  The  method  to  construct  these  fields  is 
described  in  the  next  subsection.  For  a  DNS  case  of  turbulent  pre¬ 
mixed  flame,  only  the  turbulent  velocity  field  u  is  preprocessed. 
The  species  mass  fractions  and  temperature  at  the  inflow  boundary 
are  set  to  be  constant  values  which  correspond  to  the  initial  mix¬ 
ture  composition. 

2.2.  Preprocessing  of  initial  and  inflow  mixture  fields 

Direct  numerical  simulation  of  a  complete  MILD  combustion 
system  is  not  yet  feasible  because  of  the  heavy  computational  cost 
involved.  So,  the  simulation  is  split  into  two  phases  to  mimic  the 
physical  processes  of  MILD  combustion  as  noted  in  the  previous 
subsection.  The  first  phase  involves  preprocessing  of  a  non-uni- 
form  and  inhomogeneous  mixture  field  which  is  consistent  with 
the  required  turbulence  and  combustion  conditions.  The  second 
phase  is  the  turbulent  MILD  combustion  in  which  the  inhomoge¬ 
neous  mixture  field  generated  in  the  first  phase  is  used.  These 
two  phases  are  represented  pictorially  in  Fig.  1.  Here,  the  symbol 
</>inj  is  the  equivalence  ratio  and  Q.  represents  cooling  and  heating 
of  the  product  and  reactant  mixtures  respectively.  The  desired 
mixture  fields  of  u,  7,  and  T  are  obtained  following  the  steps  in 
[12,24],  which  are  described  briefly  below. 

A  fully  developed  homogeneous  isotropic  turbulence  field  is 
obtained  first  by  simulating  a  freely  decaying  non-reacting  turbu¬ 
lence  field  evolving  from  a  specified  initial  spectrum  25].  Then,  a 
homogeneous  scalar  field  is  obtained  by  specifying  a  scalar  energy 
spectrum  as  in  [26],  which  is  taken  as  the  initial  field  of  reaction 
progress  variable  cY  defined  as  cY  =  1  -  Y//Y/,r,  varying  from  0  to 
1,  where  Yf  is  the  fuel  mass  fraction  and  the  subscript  r  denotes 
the  reactant  mixture.  The  species  mass  fractions  obtained  from  a 
one-dimensional  freely  propagating  laminar  flame  with  desired 
MILD  conditions  are  then  mapped  on  to  the  cY  field.  The  laminar 
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Fig.  1.  Schematics  of  the  two-stage  approach  employed  for  the  present  MILD  combustion  DNS.  Left  box:  mixing  DNS  domain.  Right  box:  combustion  DNS  domain.  Steps 
(i)-(iv)  noted  in  the  left  box  are  described  in  Section  2.2. 


flame  is  computed  using  the  same  chemical  kinetics  mechanism 
used  for  the  DNS.  The  unburnt  mixture  of  this  laminar  flame  is 
diluted  with  products  of  fully  burnt  mixture  (Xh2o  :  XCo2  =2:1) 
so  that  the  02  mole  fraction  in  the  reactant  mixture  X02J  matches 
a  desired  level.  The  initial  temperature  is  set  to  a  constant 
value,  Tm,  to  be  specified  later.  Finally,  these  scalar  and  velocity 
fields  are  allowed  to  evolve  for  about  one  large  eddy  turnover  time 
in  a  periodic  domain  without  any  chemical  reaction  to  mimic  EGR- 
mixing  for  a  duration  much  shorter  than  the  reference  autoignition 
delay  time,  Tign}  for  the  chosen  mixture  conditions.  The  internal 
energy  equation  is  also  solved  during  this  mixing  process,  which 
creates  a  maximum  temperature  fluctuation  of  about  2%  of  the 
mean  value,  Tm.  A  sample  field  is  shown  on  the  left  in  Fig.  1.  The 
mean  and  variance  of  the  cY  field  at  the  end  of  this  EGR  mixing  pro¬ 
cess  are  respectively  (cY)  «  0.50  and  (c? )  «  0.09  for  all  the  MILD 
cases  considered  for  this  study.  More  elaborate  discussion  on  the 
construction  of  the  initial  fields  can  be  found  in  [12,24]. 

Bilger’s  mixture  fraction  [27  has  a  variation  of  about  ±5%  of  the 
mean  value,  (£),  in  the  initial  field.  The  equivalence  ratio  obtained 
using  0  =  (1  -  £SfK/(l  -  where  £st  is  the  stoichiometric 
mixture  fraction,  gives  a  mean  value  of  (0)  =  0.8  for  all  the  cases 
considered  in  this  study.  The  calculation  of  £  is  based  on  the 
boundary  condition  for  the  air  stream  diluted  with  products  to  a 
desired  level  of  oxygen  and  a  pure  fuel  stream  as  shown  in  Fig.  1. 

A  typical  cY  field  generated  at  the  end  of  the  EGR  mixing  process 
is  shown  in  Fig.  2.  The  contours  depicted  for  a  typical  x-y  plane 
shows  that  the  field  is  not  only  inhomogeneous  but  also  turbulent. 


x 

(a) 


The  PDF  of  cY  constructed  from  the  initial  field  is  also  shown  in 
Fig.  2.  The  turbulent  mixing  and  molecular  diffusion  that  have 
occurred  during  the  EGR  mixing  process  yield  samples  with 
intermediate  values,  0  <  cY  <  1.  The  PDF  of  cY  shows  that  a 
substantial  portion  of  mixture  is  well  mixed  after  this  procedure, 
while  there  are  local  mixtures  which  still  have  either  fresh  or 
exhaust  gases.  The  PDF  is  shown  for  three  MILD  cases  considered 
in  this  study  and  the  thermo-chemical  and  turbulence  conditions 
of  these  cases  are  discussed  next. 

2.3.  Combustion  conditions 

The  thermochemical  conditions  of  the  laminar  flames  used  in 
the  generation  of  the  initial  and  inflowing  fields  described  in  Sec¬ 
tion  2.2  are  given  in  Table  1.  The  flames  A  and  B  have  diluted  reac¬ 
tant  mixture  as  signified  by  the  mole  fraction,  X,-,  r,  of  species  i  in 
the  reactant  mixture  and  thus  these  two  laminar  flames  are  used 
in  the  construction  of  mixture  fields  for  the  turbulent  MILD  cases. 
The  burning  velocity  of  these  laminar  flames  is  given  in  Table  1  as 
SL  and  the  thermal  thickness  is  defined  as  Sth  =  (Tp  -  Tr)/|VT|max. 
The  reference  autoignition  delay  time  for  these  mixtures  is  also 
given  in  the  table.  Here,  this  autoignition  delay  times  are  computed 
using  a  zero-dimensional,  constant  pressure  homogeneous  reactor 
model  and  employing  the  maximum  temperature  gradient 
criterion.  Flame  C  is  a  classical  premixed  flame,  which  is  indicated 
by  the  large  temperature  rise,  Tp  -  Tr,  and  the  reactant  mole 
fractions  listed  in  Table  1. 


Fig.  2.  (a)  Spatial  variation  of  cY  field  at  mid  x-y  plane  preprocessed  using  the  steps  explained  in  Section  2.2  for  Case  Bl.  (b)  PDF  of  cY  in  the  preprocessed  field.  Solid  line: 
preprocessed  field  to  be  used  for  Case  Al,  dashed  line:  for  Case  A2,  and  dash-dotted  line:  for  Case  Bl.  Conditions  of  these  cases  are  explained  in  Section  2.3. 
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Table  1 

Summary  of  laminar  flame  properties.  Units  are  Kelvin  (K),  metre  (m)  and  second  (s). 


XCH4.r 

X0  2,r 

^H20,r 

^C02,r 

Tr 

TP 

Sl 

103<5t„ 

103Tign 

Flame  A 

0.019 

0.048 

0.121 

0.061 

1500 

1865 

3.20 

0.69 

4.97 

Flame  B 

0.014 

0.035 

0.132 

0.066 

1500 

1775 

2.15 

0.94 

6.71 

Flame  C 

0.078 

0.194 

0.0 

0.0 

600 

2179 

1.18 

0.37 

- 

MIFE  A 

0.0095 

0.035 

0.136 

0.064 

1500 

1692 

2.62 

1.00 

5.42 

MIFE  B 

0.0060 

0.025 

0.143 

0.068 

1500 

1624 

1.66 

1.29 

6.96 

The  temperature  rise  across  the  laminar  flames  with  the  diluted 
and  preheated  mixtures  is  observed  to  be  relatively  small,  about 
370  I<  for  Flame  A  and  280  I<  for  Flame  B,  which  is  typical  of  MILD 
combustion.  Spatial  variations  of  species  mass  fraction  in  the 
inflowing  mixture  fields  are  inevitable,  as  seen  in  Fig.  2a,  when 
the  MILD  combustion  occurs  in  three  dimensions  and  it  is  not  pos¬ 
sible  to  represent  this  spatial  variation  in  a  representative  ID  lam¬ 
inar  flame.  Indeed,  the  effect  of  additional  dilution  resulting  from 
spatial  variation  reduces  the  burnt  mixture  temperature  in  the  tur¬ 
bulent  cases  by  about  150  I<  compared  to  the  respective  laminar 
flames  [12  .  To  account  for  the  dilution  effect  even  for  laminar 
flames,  the  mass  fractions  of  major  species  in  the  reactant  mixture 
are  computed  based  on  the  volume  average  of  the  respective  spe¬ 
cies  in  the  preprocessed  DNS  initial  field  obtained  after  the  EGR 
mixing  step.  A  laminar  flame  having  this  reactant  mixture  is 
named  as  a  MILD  Flame  Element  (MIFE)  in  this  study  [24  .  As 
one  would  expect,  this  averaging  yields  significantly  smaller  mole 
fraction  values  as  noted  in  Fable  1  for  MIFEs  A  and  B  compared  to 
the  Flames  A  and  B.  Also,  the  burnt  mixture  temperature,  Tp,  is  only 
about  27  K  higher  than  the  DNS  values  (this  is  the  maximum  differ¬ 
ence  observed).  For  these  reasons,  MIFEs  are  taken  to  be  represen¬ 
tative  flamelets  for  the  MILD  combustion  cases  investigated  in  this 
study  and  their  thermochemical  attributes  listed  in  Table  1  are 
used  to  normalise  the  respective  DNS  results.  For  example,  length, 
gradient  of  reaction  progress  variable  and  reaction  rate  are  respec¬ 
tively  normalised  using  Sth,!/Sth  and  prSL/Sth.  The  normalised 
quantities  are  denoted  using  a  superscript  “+”. 


23  A.  Conditions  for  turbulent  MILD  combustion 

Three  MILD  and  one  classical  premixed  combustion  case  were 
simulated.  The  turbulence  and  thermochemical  conditions  of  these 
four  cases  are  given  in  "able  2.  The  thermochemical  conditions  for 
Cases  Al  and  A2  are  based  on  Flame  A.  Thus,  they  have  the  same 
dilution  level,  but  their  turbulence  conditions  are  different  as 
shown  in  "able  2.  Case  Al  has  the  same  turbulence  field  as  in  Case 
Bl,  but  uses  the  same  thermochemical  conditions  as  Case  A2.  All  of 
the  four  cases  have  an  equivalence  ratio  of  </>  =  0.8  and  the  autoig¬ 
nition  temperature  for  this  lean  methane-air  mixture  is  1100K. 
The  temperature  of  the  initial  and  inflowing  mixture  is  set  to  be 
Tm  «  1500  I<  for  the  MILD  combustion  cases,  and  600  I<  for  the 
premixed  combustion  Case  C.  The  high  value  of  Tm  for  the  MILD 
combustion  cases  is  comparable  to  that  used  in  experiments 
[28].  This  inlet  temperature  together  with  the  low  mole  fraction 
of  02  suggests  that  the  combustion  conditions  are  strictly  in  the 
MILD  regime. 


The  maximum  and  averaged  (Xo2,r)  mole  fractions  of  oxy¬ 
gen  in  the  reactant  mixture  for  the  MILD  cases  are  significantly 
smaller  than  for  Case  C  because  of  the  high  dilution  levels  (see 
Table  2).  These  dilution  levels  are  comparable  to  previous  studies 

[2,6,3]. 

The  values  of  mean  (£)  and  stoichiometric  £st  mixture  fractions 
are  also  given  in  fable  2.  The  RMS  of  velocity  fluctuations  and  the 
integral  length  scale  of  the  initial  turbulence  field  are  denoted 
respectively  as  u!  and  l0  in  Table  2.  The  Zeldovich  thickness  is 
dF  =  (2/pcp)/SL,  where  2  and  cp  are  respectively  the  mixture  ther¬ 
mal  conductivity  and  specific  heat  capacity  at  constant  pressure. 
The  Reynolds  number,  Re/o  =  u'l0/vr ,  with  vr  as  the  kinematic  vis¬ 
cosity  of  reactant  mixture,  is  varied  respectively  from  67  to  96 
for  the  MILD  cases.  An  estimate  of  the  Reynolds  number  for  MILD 
combustion  in  previous  studies  [29,11,13]  using  the  information 
for  reacting  jet  flows  [30-33]  suggests  that  the  values  of  Re/o  in 
Table  2  are  comparable  to  those  for  previous  MILD  combustion 
experiments.  The  turbulence  level  for  Case  C,  the  premixed 
combustion  case,  is  deliberately  set  to  be  a  small  value  to  help  to 
contrast  the  behaviour  of  reaction  zones  and  their  structure  in  MILD 
conditions  and  in  premixed  combustion  with  large  Damkohler 
number.  The  Damkohler  and  Karlovitz  numbers  are  defined  as 
Da  =  ( lo/SF)/(w/SL )  and  Ka  «  (u//SL)3/2(/0/<5F)“1/2.  The  three  MILD 
cases  are  in  the  thin-reaction  zones  regime,  and  Case  C  is  near  the 
border  between  the  thin-reaction  zones  and  corrugated  flamelets 
regimes  in  a  turbulent  combustion  regime  diagram  [34]. 

2.3.2.  Computational  detail 

The  computational  domain  is  of  size  Lx  x  Ly  x  Lz  =  lO.Ox 
10.0  x  10.0  mm3  for  the  three  MILD  cases  and  Lx  x  Lyx 
Lz  =  10.0  x  5.0  x  5.0  mm3  for  the  premixed  flame,  Case  C.  If  the 
computational  domain  is  normalised  using  3th  of  the  respective 
laminar  flames  (MIFEs  for  MILD  and  Flame  C  for  premixed  cases), 
If  x  Ly  x  If  =  103  for  Cases  Al  and  A2,  (7.75)3  for  Case  Bl,  and 
26.7  x  13.4  x  13.4  for  Case  C.  These  domains  are  discretised 
using  512  x  512  x  512  mesh  points  for  Cases  Al  and  A2, 
384  x  384  x  384  mesh  points  for  Case  Bl,  and  512  x  256  x  256 
mesh  points  for  Case  C.  These  meshes  ensure  that  there  are  at  least 
15  mesh  points  inside  the  respective  thermal  thickness  Sm. 

The  simulations  of  the  MILD  cases  were  run  for  1.5 
flow-through  times  to  ensure  that  the  initial  transients  had  left 
the  domain  before  collecting  data  for  statistical  analysis.  The 
flow-through  time  is  defined  as  td  =  Lx/Uin,  which  is  the  mean 
convection  time  from  the  inflow  to  the  outflow  boundary.  The 
simulations  were  then  continued  for  one  additional  flow-through 
time  and  80  data  sets  were  collected.  For  Case  C,  93  data  sets  were 


Table  2 

Turbulent  combustion  conditions. 


vmax 

A02,r 

<*02,r> 

(0 

Zst 

Uin/SL 

u'/Sl 

lo/$F 

lo /$th 

Re/o 

Da 

Ka 

Case  Al 

0.048 

0.035 

0.011 

0.014 

9.6 

6.26 

10.8 

1.48 

96.2 

1.72 

4.78 

Case  A2 

0.048 

0.035 

0.011 

0.014 

9.6 

3.80 

12.3 

1.70 

67.0 

3.25 

2.11 

Case  Bl 

0.035 

0.025 

0.008 

0.010 

15.1 

9.88 

6.8 

1.15 

96.1 

0.69 

11.9 

Case  C 

0.194 

0.194 

0.044 

0.055 

3.0 

2.19 

12.3 

2.11 

38.5 

5.64 

0.92 
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collected  over  a  time  of  0.56td  after  allowing  one  flow-through 
time  for  initial  transients  to  exit  the  computational  domain.  These 
simulations  have  been  run  on  Cray  XE6  systems  using  4096  cores 
with  a  wall-clock  time  of  about  120  h  for  Case  Al  and  Case  A2, 
which  have  the  largest  number  of  mesh  points,  and  using  16,384 
cores  with  80  h  of  wall-clock  time  for  Case  C. 

3.  Results  and  discussion 

3.1.  Temperature  and  reaction  zone  behaviour 

Typical  PDFs  of  normalised  temperature,  cT  =  (T  -  Tr)/(TP  -  Tr), 
are  shown  in  Fig.  3  for  various  x+  locations  for  the  Case  B1  in 
Table  2.  These  PDFs  are  constructed  using  samples  collected  over 
the  entire  sampling  period.  This  PDF  indicating  a  bimodal  behav¬ 
iour  for  the  premixed  case  is  also  shown  for  x+  =  17.5  in  Fig.  3. 
The  PDFs  of  the  MILD  case  show  neither  a  sharp  peak  nor  a  bimodal 
behaviour.  The  PDF  for  the  x+  =  0.30  location,  which  is  near  the 
inlet  boundary,  shows  that  unburnt  gases  are  predominant  at  this 
region.  The  PDF  shows  a  plateau  in  0.2  ^  cT  <  0.7  for  the  location 
x+  =  2.9  suggesting  distributed  combustion.  As  one  moves  further 
downstream,  the  probability  of  finding  burnt  gases  increases  and 
the  PDF  distribution  is  relatively  broader  compared  to  that  of  the 
premixed  case. 

Typical  contours  of  cT  and  cbCr,  which  is  the  normalised  reaction 
rate  of  cT,  are  shown  in  Fig.  4  for  a  x-y  plane  located  at  z  =  Lz/2.  The 
reaction  rate  is  calculated  as  cdCt  =  Q/cp(Tp  -  Tr),  where  Q  is  the 
heat  release  rate.  The  contours  of  coCt  =  coCT/coCT,max  =  0.2  and  0.5 
shown  in  Fig.  3  for  the  premixed  case  form  thin  reaction  zones  with 
a  typical  thickness  of  about  dth,  where  coCr,max  is  the  global  maxi¬ 
mum  reaction  rate  at  t  =  1.5td.  This  behaviour  is  consistent  with 
the  bimodal  PDF  shown  in  Fig.  3  for  the  premixed  case.  Further¬ 
more,  the  peak  reaction  rate  is  located  along  cT  =  0.6  contour  (grey 
dashed  line),  which  is  consistent  with  a  previous  observation  [35] 
on  the  relation  between  cT  and  coCt.  Although  the  value  of  Da  for 
Case  C  is  larger  than  for  the  MILD  cases  (see  Table  2),  previous 
DNS  [36-38  of  turbulent  premixed  flames  with  Da  <  0.69  showed 
thin  sheet-like  reaction  zones  with  some  interactions. 

For  the  MILD  cases,  regions  with  coCt  ^  0.2  occupy  a  substantial 
part  of  the  combustion  domain  and  the  extent  of  such  regions 
increases  with  dilution  and  turbulence  levels.  Thin  reaction  zones 
are  also  observed  locally,  which  are  typically  surrounded  by  the 
C0cT  ^  0.5  contours.  Such  regions  do  not  seem  to  follow  a  particular 
cT  contour  as  for  the  premixed  case.  Also  these  regions  are  spread 
out  in  the  domain  to  form  a  patchy  appearance.  Such  patchy 
reaction  zones  are  also  observed  in  OH  PLIF  images  of  a  previous 


Fig.  3.  PDFs  of  cT  for  various  x+  locations  for  Case  Bl.  A  typical  PDF  of  the  premixed 
case  is  also  shown  for  x+  =  17.5. 


experimental  study  [8  .  Despite  the  existence  of  such  thin  local 
reaction  zones,  the  temperature  field  is  distributed  as  shown  in 
the  PDF  of  cT  in  Fig.  3.  This  results  from  broader  and  distributed 
reaction  zones,  which  are  also  seen  in  Fig.  4.  The  patchy  appear¬ 
ance  or  distributed  nature  of  the  reaction  zones  results  from  their 
abundant  and  frequent  interactions.  Previous  analysis  of  time  his¬ 
tories  of  the  reaction  zones  has  identified  these  interactions  [12], 
which  are  not  considered  further  in  this  study.  Also,  the  intense 
reaction  zones  near  the  inlet  boundary  observed  in  Fig.  4  for  the 
MILD  cases  are  due  to  the  radicals  and  intermediate  species  pres¬ 
ent  in  the  reactant  mixture  diluted  using  products.  The  presence 
of  these  species  and  the  high  temperature  of  the  reactant  mixture 
will  significantly  shorten  the  autoignition  delay  time. 

The  reaction  zones  represented  by  cbCr  ^  0.5  in  Fig.  4  are  ana¬ 
lysed  next  to  study  their  morphological  characteristics  in  order 
to  understand  whether  these  reaction  zones  are  thin  or  not. 


3.2.  Morphology  of  reaction  zones 


3.2.1.  Minkowski  functional  and  shapefinders 

The  reaction  zone  morphology  is  analysed  using  Minkowski  func¬ 
tionals  [39],  which  help  to  define  shapefinders.  The  shapefinders  are 
characteristic  scales  for  length,  L,  width,  W,  and  thickness,  T,  of  a 
three  dimensional  object  of  interest.  The  Minkowski  functionals 
and  shapefinders  have  been  used  to  quantify  iso-density  structures 
in  cosmology  [40-48],  magnetic  field  structures  in  magneto¬ 
hydrodynamics  [49]  and  structures  of  homogeneous  turbulence 
[50].  There  are  (n  +  1)  Minkowski  functionals  for  a  n-dimensional 
object  and  the  4  functionals  for  a  3D  object  are  given  as  [43,48] 


Vo 

Vi 

V2 

V3 


V, 

s 

6’ 

1 

3  n 

1 

2n 


K\  +  K2 
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(k\K2)  ds , 


a) 

(2) 

(3) 

(4) 


where  V  is  the  volume  of  the  object  enclosed  by  its  surface,  S,  with  a 
surface  area  of  S ,  and  K\  and  k2[k\  ^  k2)  are  the  two  principal  cur¬ 
vatures  at  a  given  point  on  S.  It  is  apparent  that  the  mean  and 
Gaussian  curvatures  integrated  over  the  entire  surface  are  involved 
respectively  in  V2  and  \Z3,  and  local  values  of  these  curvatures  are 
used  in  many  turbulent  combustion  models.  The  integrated 
Gaussian  curvature,  V3,  is  the  Euler  characteristic,  y,  which  is  a 
topological  invariant  of  the  three-dimensional  object  of  interest. 
These  four  functionals  defined  in  Eqs.  ( 1 )— (4)  are  Galilean  invariant 
morphological  properties  of  a  three  dimensional  object.  In  the  pres¬ 
ent  study,  the  Minkowski  functionals  are  computed  for  many 
objects  identified  using  iso-surfaces  of  reaction  rate  and  extracted 
by  employing  a  grid  cell  counting  method  [42,44]. 

Three  characteristic  scales  for  the  shapefinder  are  obtained 
from  the  four  functionals  given  in  Eqs.  (l)-(4)  as  [43]: 


Thickness,  T  =  ^ - , 

2Vi 

(5) 

Width,  W  ~  2  , 
nV2 

(6) 

3  VG 

Length,  L  =  2V3, 

(7) 

and  they  satisfy  L  ^  W  ^  T  for  any  convex  body  [44].  Also  this 
inequality  holds  for  closed  partially  concave  objects  since  the 
Minkowski  functionals  are  additive  [44,50  .  These  four  functionals 
have  to  be  positive  to  obtain  physically  meaningful  shapefinders 
using  Eqs.  (5)-(7).  However,  V3  =  %  can  be  positive  or  negative  or 
zero  because  it  is  related  to  the  three-dimensional  genus 
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(a) 


(b) 


(c) 


(d) 


Fig.  4.  Contours  of  cT  and  cjCt  in  the  mid  x-y  planes  for  (a)  Case  C,  (b)  Case  Al,  (c)  Case  A2,  and  (d)  Case  Bl.  Thin  solid  line  is  coCt  =  0.2,  thick  solid  line  is  coCt  =  0.5,  and  thick 
grey  dashed  line  is  cT  =  0.6.  Note  that  only  5.3  <  x+  ^  22.2  of  the  computational  domain  is  shown  for  Case  C.  All  snapshots  are  taken  at  t  =  1.5td. 


G  =  1  -  0.5/,  which  is  defined  as  the  number  of  cuts  that  can  be 
made  along  a  simple  curve  on  an  object  without  splitting  it.  For 
an  object  made  by  a  single  closed  surface  G  =  0  and  thus  V3  =  2. 
If  there  is  a  hole  on  the  object,  say  for  example  a  torus,  then 
G  =  1  and  thus  V3  =  0.  For  a  pretzel  G  =  2  and  thus  V3  =  -2.  Hence, 
V3  is  related  to  the  number  of  holes,  Nhoie,  on  the  object  as 
V3  w  1  -  Nh()ie.  This  creates  a  problem  for  the  definition  of  L  in  Eq. 
(7)  because  L  is  singular  for  V3  =  0  and  it  becomes  unphysical  when 
V3  <  0.  To  avoid  this,  an  alternative  definition  of  L  is  written  for 
objects  with  non-zero  G  or  V3  ^  0  as  [45]: 


I 


3V2 

4(gTT)- 


The  definition  of  L  in  Eqs.  (7)  and  (8)  is  identical  for  objects  with  no 
holes,  G  =  0.  For  objects  with  holes,  G  >  0,L  signifies  a  characteristic 
length  between  holes  and  it  implies  a  representative  radius  of  a  cir¬ 
cular  torus  when  G  =  1  [50].  For  a  sphere  of  radius  r,  Eqs.  (l)-(8) 
give  T  =  W  =  L  =  r.  For  an  infinitely  long  cylinder  with  radius  X\ , 
these  length  scales  are  (T,  W,L)  =  (1.5rl52ri,oo)  and  for  an  infi¬ 
nitely  large  thin  circular  disk  of  thickness  3,  they  are 
(T,  W,  L)  =  (1.5 3  ,  oc,  oo).  Thus,  these  three  quantities  yield  represen¬ 
tative  scales  for  the  spatial  extent  of  a  structure  and  do  not  give  its 
exact  dimensions  except  for  a  sphere. 

Dimensionless  shapefinders  are  defined  as  [43]: 

W  -T 

Planarity,  P  =  ,  (9) 

Filamentarity,  F  =  j-  -  ^ .  (10) 


These  two  quantities  allow  a  convenient  visual  representation  of 
the  geometry  for  a  group  of  objects  in  a  two-dimensional  (P,  F) 
space.  They  vary  from  0  to  1  when  the  inequality  L  ^  W  ^  T  is 
met.  The  filamentarity  approaches  to  1  when  L  is  very  large  com¬ 
pared  to  W  and  T.  The  planarity  is  close  to  1  when  W  >  T.  Thus  a 
long  tube  has  F  ~  1  and  P  ~  0  while  a  sheet  has  F  ~  0  and  P  ~  1. 
Typical  shapes  are  marked  in  the  shapefinder,  (P,F),  map  shown 
in  Fig.  5.  It  is  obvious  that  (P,  F)  =  (0, 0)  for  a  spherical  object. 


3.2.2.  Shape  of  reaction  zones 

The  reaction  zones  are  identified  first  using  a  threshold  for  co+ 
and  this  identification  introduces  some  subjectivity  in  the  analysis 

Infinitely  long  Infinitely  long 


Fig.  5.  A  representative  map  of  shapefinders  [50]. 
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depending  on  the  threshold  value  used.  Also,  the  use  of  a  single 
threshold  value  for  all  the  cases  does  not  identify  regions  with  sim¬ 
ilar  reaction  rate  magnitudes  since  the  normalising  factor  prsL/Sth 
is  not  the  same  for  all  the  cases  considered  in  this  study.  To  remove 
these  ambiguities,  the  threshold  value  is  chosen  using  Z(g)  = 
(A SNC  cOct\cdCj  =  C)  rather  than  just  cdCt,  where  Nc  =  Dc(Vcr  •  VcT) 
is  the  scalar  dissipation  rate  and  AS  is  the  surface  area  identified 
using  coCt  =  The  scalar  dissipation  rate  and  reaction  rate  are  large 
when  the  reaction  zones  have  propagative  characteristics  (rather 
they  are  propagating  flames).  Thus,  a  large  value  of  Z  automati¬ 
cally  implies  propagating  flames  and  so  their  morphological 
characteristics  can  be  compared  meaningfully  to  those  for 
premixed  flames.  Furthermore,  using  only  a  threshold  value  for 
does  not  allow  us  to  differentiate  between  the  reaction  and 

Ct 

propagating-flame  dominated  structures. 

The  variations  of  Z(£)/Zmax  with  £+  for  all  the  cases  are  shown 
in  Fig.  6  and  the  value  of  £+  corresponding  to  Z(£)/Zm ax  =  1  in  each 
case  is  chosen  as  the  threshold  to  identify  reaction  zone  iso-sur- 
faces  for  morphological  analysis.  It  is  found  that  these  surfaces 
enclose  about  the  top  50%  of  the  reaction  rate  in  each  case  and  thus 
they  represent  regions  with  intense  reaction.  It  has  also  been  ver¬ 
ified  that  the  statistics  of  shapefinders  do  not  change  unduly  when 
the  threshold  is  set  to  a  value  20%  smaller  than  fj\  Furthermore, 
the  use  of  cY  and  its  reaction  rate  does  not  change  the  conclusions 
from  this  analysis.  There  were  unclosed  objects  because  of  their 
interaction  with  inflow/outflow  boundaries  and  these  objects  were 
excluded  from  further  consideration  since  the  Minskowski  func¬ 
tionals  require  closed  objects. 

For  most  of  the  objects  identified  using  the  threshold  value  of  £+ 
as  noted  above,  the  surface  average  of  K\  k2  is  positive  for  all  the 
four  cases  considered  in  this  study.  This  implies  that  V3  is  also 
positive  for  these  objects.  However,  the  joint  PDFs  of  K\  and  k2 
(not  shown  here)  for  the  identified  structures  suggest  that  there 
are  some  negative  principal  curvature  (s)  locally.  Thus,  these  struc¬ 
tures  have  local  concavity  possibly  with  holes,  although  they  are 
made  of  convex  surface  in  an  integral  sense.  Therefore,  Eq.  (8)  is 
used  to  compute  L  in  the  present  analysis.  Also,  a  small  fraction 
(about  15%)  of  these  structures  have  a  multiply  connected  surface 
because  of  the  presence  of  holes  resulting  from  reaction  zone  con¬ 
tortions  and  their  interaction  due  to  turbulence. 

The  objects  extracted  using  the  post  processing  steps  discussed 
above  are  analysed  further  to  understand  if  the  MILD  combustion 
involves  flame  sheets,  thin  or  distributed  reaction  zones.  Two  of 
these  objects  chosen  arbitrarily  are  shown  in  Fig.  7  as  examples 
and  these  are  obtained  from  the  MILD  cases  Al,  Fig.  7a,  and  A2, 
Fig.  7b.  As  one  sees,  these  objects  are  made  of  multiply  connected 
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Fig.  6.  Variation  of  Z(£)/Zmax  with  £+  in  Case  Al  (solid  line),  Case  A2  (dashed),  Case 
B1  (dash-dotted)  and  Case  C  (dotted). 


surface  and  it  is  not  easy  to  say  if  they  can  be  classified  as  thin 
sheets  or  thin  reaction  zones  based  on  a  visual  examination  with¬ 
out  any  subjectivity.  The  Minkowski  functionals  and  shapefinders 
discussed  in  the  previous  subsection  help  us  on  this  matter  as 
one  shall  see  in  the  following  discussion. 

The  values  of  these  functionals,  three  characteristic  length 
scales  and  shapefinders  are  given  in  the  caption  of  Fig.  7,  and  they 
satisfy  the  inequality  L  ^  W  ^  T.  The  object  shown  in  Fig.  7a  pre¬ 
sents  five  holes,  indicated  by  the  arrows,  and  this  implies  that 
V3  <  0  as  noted  in  the  previous  subsection.  In  this  case,  Eq.  (8)  is 
used  to  compute  L.  The  values  of  shapefinders,  P  and  F,  clearly  sug¬ 
gests  that  these  objects  are  pancake-like  in  the  shapefinders  map 
shown  in  Fig.  5.  Thus,  the  morphology  of  the  complex  shaped 
reaction  zones  can  be  identified  unambiguously  and  analysed 
statistically. 

Figure  8  shows  arbitrarily  chosen  typical  shapes  of  reaction 
zones  along  with  the  values  of  their  shapefinders,  (P,F),  for  the 
MILD  and  premixed  cases  considered  here.  This  analysis  is  con¬ 
ducted  using  half  the  number  of  data  sets  collected  and  including 
the  full  data  sets  does  not  changes  the  statistics  to  be  discussed 
in  the  remainder  of  this  subsection  unduly.  There  are  510,  328, 
272  and  42  objects  respectively  for  Cases  Al,  A2,  B1  and  C  for 
the  statistical  analysis  of  morphological  characteristics.  Since  the 
reaction  rate  threshold,  ,  used  is  relatively  large,  even  the  statis¬ 
tically  planar  premixed  flame  has  objects  enclosed  by  surfaces  not 
interacting  with  boundaries.  These  objects  are  sheet-like  as  in 
Fig.  8f-j  for  the  premixed  case  and  their  shapes  are  varied,  complex 
and  multiply  connected  for  the  MILD  combustion  cases.  It  is  worth 
noting  that  there  are  spherical  reaction  zones  in  the  MILD  combus¬ 
tion  as  in  Fig.  8c  clearly  indicated  by  the  low  (P,F)  values. 

Figure  9  shows  the  three  length  scales,  L,  W  and  T,  normalised 
using  the  respective  Sth  for  the  MILD  and  premixed  cases.  The  black 
solid  lines  are  for  W+  =  T+  and  L+  =  W+  and  are  used  to  verify 
whether  the  inequality  T+  ^  W+  ^  L+  is  met  by  the  computed 
results  in  order  to  have  meaningful  P  and  F.  As  one  observes  in 
Fig.  9,  this  inequality  is  satisfied  predominantly  and  there  are  some 
samples  violating  the  condition  L+  >  W+.  A  detailed  investigation 
of  these  samples  identified  that  this  is  due  to  numerical  errors  as 
these  objects  are  small  and  they  are  excluded  from  further  consid¬ 
eration  since  they  give  negative  F  (see  Eq.  (10))  which  is  non-phys¬ 
ical.  The  effect  of  spatial  resolution  on  the  calculation  of  the 
shapefinders  is  explored  in  Appendix  A.  Furthermore,  it  has  been 
observed  that  only  a  very  small  percentage  of  the  reaction  zone 
structures  extracted  from  the  DNS  does  not  obey  the  inequality. 
These  percentages  are  1.8%  in  Case  Al,  1.2%  in  Case  A2,  0.37%  in 
Case  B1  and  none  in  Case  C.  Thus,  it  becomes  clear  that  the  statis¬ 
tics  reported  here  are  not  influenced  by  omitting  these  small  struc¬ 
tures  or  numerical  errors  while  calculating  P  and  F. 

Figure  9  shows  that  W+  is  much  larger  than  T+  and  it  is  compa¬ 
rable  to  L+  for  Case  C,  which  is  the  premixed  combustion  case.  This 
variation  suggests  that  the  extracted  reaction  zone  structures  are 
thin  sheets  as  one  would  expect  in  premixed  combustion.  The 
characteristic  thickness  for  the  premixed  reaction  zones  vary  from 
0.05  ^  T+  ^  0.15,  which  is  relatively  narrow.  Their  characteristic 
length  and  width,  however,  vary  over  a  wide  range  as  shown  in 
Fig.  9b  and  the  data  follow  W  =  L  line  closely,  suggesting  a  thin 
sheet.  The  range  for  these  characteristic  scales  is  wider  for  the 
MILD  combustion  cases  as  one  sees  in  Fig.  9.  The  typical  ranges 
are  0.09  ^  F+ ^  0.5,  0.2  ^  W+  ^  1.2  and  0.2  ^  L+  ^  5.5.  These 
ranges  are  reduced  relatively  in  Cases  A2  and  B1  compared  to  Case 
Al.  As  noted  in  "able  2,  Cases  A2  has  a  lower  turbulence  Reynolds 
number  compared  to  Case  Al  despite  the  same  dilution  level  and 
thus  the  range  of  these  scales  are  reduced  somewhat.  However, 
for  the  most  diluted  Case  B1  having  the  same  Rei0  as  Case  Al,  most 
of  the  data  are  concentrated  in  a  relatively  narrower  range.  This  is 
because  of  mixture  non-uniformity  in  Case  B1  resulting  in 
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(a)  (b) 


Fig.  7.  Two  objects  chosen  arbitrarily  are  shown  above  and  arrows  indicate  the  locations  of  holes  on  them.  The  values  of  their  Minkowski  functionals  and  shapefinders  are  as 
follows.  For  the  extracted  object  in  (a)  V0  =  8.1  mm3,  Vi  =  14.5  mm2,  V2  =  11.0  mm,V3  =  -8.0  (G  =  5 ),T  =  0.28  mm,  W  =  0.84  mm,L  =  1.4  mm,  (P,F)  =  (0.50,0.25).  For 
the  object  in  (b)  V0  =  4.6  mm3,  V\  =  9.1  mm2,  V2  =  6.7  mm,l/3  =  0.0  (G  =  !),!  =  0.26  mm,  W  =  0.87  mm,  l  =  2.5  mm,  (P, F)  =  (0.54,0.49). 


(b)  (0.35,0.69). 


(c)  (0.05,0.05). 


(d)  (0.48,0.42). 


(e)  (0.46,0.62). 


(a)  (0.44,0.15). 


(f)  (0.68,0.24). 


(g)  (0.70,0.18). 


(h)  (0.70,0.19). 


(i)  (0.81,0.19). 


CD  (0.67,0.23). 


Fig.  8.  Samples  of  reaction  zone  shapes  extracted  from  the  DNS  of  Case  B1  (a  to  e)  and  Case  C  (f  to  j)  at  an  arbitrary  time.  The  values  of  their  shapefinders,  (P,  F),  are  also  given 
above. 


Fig.  9.  The  three  length  scales,  calculated  using  Eqs.  (5),  (6),  and  (8),  for  the  reaction  zones  in  Case  Al  (o),  Case  A2  (A),  Case  B1  (□)  and  Case  C  (x). 


convolutions  and  interactions  of  reaction  zones  as  shown  in  Fig.  4. 
The  characteristics  of  sheet-like  structures  seen  in  the  premixed 
case  is  not  observed  for  the  MILD  combustion  cases. 

The  characteristic  scales  shown  in  Fig.  9  are  used  to  calculate 
the  shapefinders  P  and  F  using  Eqs.  (9)  and  (10),  which  are  shown 
in  Fig.  10.  The  results  for  the  premixed  case  reside  in  a  narrow 


region  mostly  confined  to  the  (P,F)  regions  for  sheet-like  struc¬ 
tures,  but  the  MILD  reaction  zones  have  a  broader  distribution  with 
about  0.01  ^  P  ^  0.7  and  0.01  ^  F  ^  0.8.  A  comparison  of  these 
results  to  Fig.  5  suggests  that  the  MILD  reaction  zone  shapes  range 
from  blobs  (small  P  and  F)  to  pancakes  (intermediate  P  and  F). 
There  are  some  short  ribbons  also  and  the  pocket  burnouts  would 
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Fig.  10.  Shapefinder,  P  and  F,  for  Case  Al  (o),  Case  A2  (A),  Case  B1  (□),  and 
Case  C  (x). 


larger  volumes  in  Case  Al  due  to  frequent  reaction  zone  interac¬ 
tions  resulting  from  higher  turbulence  level  although  these  two 
cases  have  equal  number  (in  a  probability  sense)  patchy  reaction 
zones  having  V+  ~  0.07. 

Generally  L+ean  is  larger  than  W+ean  and  the  ratio  L^em/W^em  is 
about  2.3  for  the  MILD  cases  and  1.5  for  the  premixed  case.  The 
ratio  W^ean/T^ea n  is  also  about  2.1  for  the  MILD  cases,  but  it  is 
7.6  for  the  premixed  case.  These  values  clearly  agree  with  the  view 
that  the  MILD  reaction  zones  are  not  thin  sheets  even  in  an  average 
sense.  Although  the  local  reaction  zone  volume  V+  would  change 
with  the  threshold  value  the  computed  value  of  V+  varies  due 
to  dilution  level  under  the  same  turbulence  condition  (compare 
Cases  Al  and  Bl).  This  suggests  that  the  widely  used  Eddy  Dissipa¬ 
tion  Concept  (EDC)  model  could  be  improved  further  for  MILD 
combustion  if  the  volume  of  local  reacting  regions  is  determined 
not  only  by  the  size  of  turbulence  energy  dissipation  regions,  but 
also  by  considering  the  effect  of  dilution  on  characteristic  scales 
such  as  V+. 


yield  a  blob-like  shape  represented  by  the  points  close  to  (0,  0)  in 
Fig.  10.  Representative  pictures  of  these  structures  as  extracted 
from  the  DNS  data  are  shown  in  Fig.  8. 

The  marginal  PDF  of  P  and  F  constructed  using  20  bins  for  the 
MILD  cases  and  12  bins  for  the  premixed  case  shown  in  Fig.  1 1  sug¬ 
gests  that  P  and  F  are  distributed  from  about  0.01  to  0.8  for  the 
MILD  cases  as  discussed  above.  However  for  the  premixed  case, 
their  ranges  are  0.6  ^  P  ^  0.9  and  0.1  ^  F  ^  0.4,  with  pronounced 
peaks  in  the  PDFs  as  seen  in  Fig.  11.  This  figure  suggests  that  the 
most  probable  value  for  (P,F)  is  (0.79,  0.2)  for  the  premixed  case. 
The  most  probable  value  of  P  varies  between  0.3  and  0.5,  and  F  is 
about  0.2  for  the  MILD  cases.  These  values  suggest  that  the  most 
probable  shape  for  the  MILD  reaction  zone  is  pancake-like.  Such 
shapes  can  result  from  autoigniting  regions  which  do  not  necessar¬ 
ily  result  in  thin  flame  fronts.  The  average  and  most  probable  val¬ 
ues  of  volume  V  and  the  three  characteristic  length  scales  of  these 
structures  used  in  the  above  analyses  are  given  in  Table  3.  These 
values  are  normalised  using  the  respective  Sth  given  in  the  same 
table.  The  mean  and  most  probable  T+  of  the  reaction  zones  show 
a  very  weak  sensitivity  to  Da,  meaning  the  dimensional  thickness 
increases  with  dilution  level  since  Sth  increases.  The  values  of 
W+ean  and  L+ean  also  show  this  insensitivity  in  the  MILD  cases. 
The  value  of  V+ean  for  Case  Al  is  about  twice  as  large  as  for  Case 
A2,  while  the  values  of  Vpro5  are  identical  for  these  two  cases.  This 
implies  that  a  substantial  number  of  local  reaction  zones  have 


3.3.  Autoignition  and  flame  propagation  characteristics 

The  analysis  in  the  previous  subsection  showed  that  MILD  reac¬ 
tion  zones  have  statistically  pancake-like  shapes.  Since  propagat¬ 
ing  flame  is  usually  associated  with  thin  sheet-like  structures, 
these  pancake-like  shapes  are  suggestive  of  frequent  events  of 
autoignition  of  interacting  flames  in  MILD  combustion.  Autoigni¬ 
tion  is  sensitive  to  the  local  temperature  and  mixture  composition 
which  are  influenced  by  turbulent  mixing  and  molecular  diffusion 
[3,51,17].  Especially,  turbulent  mixing  directly  affects  the  autoigni¬ 
tion  delay  time  [51  .  In  the  present  configuration,  inlet  mixture  is 
imperfectly  mixed  containing  reactants  and  products  of  EGR  mix¬ 
tures  and  turbulence  resulting  in  spatio-temporal  variation  of  local 
mass  fractions  of  chemically  active  radicals  and  intermediates.  This 
results  in  spatial  variation  of  autoignition  delay  time.  However,  the 
local  turbulence  conditions  influence  the  local  residence  time  for 
the  mixture  which  may  result  in  insufficient  time  for  autoignition 
but  may  create  favourable  conditions  for  propagating  flames  or 
vice-versa.  The  influences  of  global  mixture  composition  on  autoig¬ 
nition  have  been  investigated  for  MILD  combustion  in  the  past 
[52].  However,  competing  effects  of  local  turbulence  and  mixture 
composition  on  autoignition  or  propagating  flame  characteristics 
have  not  been  studied  although  these  characteristics  are  likely  to 
determine  the  type  of  modelling  used  for  MILD  combustion.  The 
predominant  characteristics  of  local  combustion  in  MILD  mixtures 
is  studied  using  a  flux-balance  analysis  in  the  next  subsection. 


Fig.  11.  PDF  of  (a)  planarity  and  (b)  filamentarity  for  Case  Al  (thin-solid),  Case  A2  (dashed),  Case  Bl  (dash-dotted),  and  Case  C  (thick  solid). 
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Table  3 

Average  and  most  probable  values  of  V  and  the  three  characteristic  length  scales  (Eqs.  (5)-(8)). 


8th  (mm) 

Da 

Ka 

V+ 

y  mean 

T+ 

1  mean 

W+ 

v  v  mean 

l+ 

Lmean 

^  pro  Id 

T+ 

1  prob 

Wjrob 

j  + 

^prob 

Case  Al 

1.00 

1.72 

4.78 

0.93 

0.19 

0.44 

1.00 

0.07 

0.15 

0.26 

0.46 

Case  A2 

1.00 

3.25 

2.11 

0.42 

0.19 

0.38 

0.87 

0.07 

0.18 

0.29 

0.45 

Case  Bl 

1.29 

0.69 

11.9 

0.30 

0.17 

0.36 

0.85 

0.12 

0.16 

0.36 

0.63 

Case  C 

0.37 

5.64 

0.92 

0.56 

0.11 

0.84 

1.27 

0.22 

0.15 

0.83 

1.23 

3.3. 1 .  Balance  of  convection,  diffusion  and  chemical  reaction 

Turbulent  combustion  involves  a  balance  of  convection,  diffu¬ 
sion  and  chemical  reaction  processes.  The  transport  equation  for 
normalised  temperature,  cT,  is  written  as: 

pDCrg)+  ,  (11) 

— v - y  U:  reaction 

C:  convection  v,  diffusion 

where  u,  is  the  fluid  velocity  in  direction  x,  and  DCr  is  the  molecular 
diffusivity  of  cT.  The  balance  among  these  terms,  C,  V  and  7 Z,  may  be 
assessed  by  studying  B  defined  as: 

B=\C-V\-\7Z\:  (12) 

which  makes  Eq.  (11)  read  as  dp  cT/dt  +  B  =  0.  For  zero-dimen¬ 
sional  homogeneous  transient  combustion  in  a  fixed-mass  constant 
pressure  reactor,  Eq.  (11)  is  dominated  by  the  chemical  source  1Z 
and  unsteady  terms,  giving  B  =  -\R\.  This  situation  is  akin  to  com¬ 
bustion  in  a  Perfectly  Stirred  Reactor  (PSR)  with  a  volume  V  and  a 
mass  flow  rate  of  rhin,  having  a  governing  equation  for  cT  given  by 

dcj  _  Cj  coCj 

~dt~~TR  p ’  1  j 

where  tr  is  the  residence  time  of  a  reactor  defined  as  tr  =  pV/min. 
Thus,  for  a  model  configuration  such  as  homogeneous  reactor  or 
PSR,  regions  with  B  <  0  signify  autoigniting  events,  where  the 
transport  equation  of  cT  is  dominated  by  the  unsteady  and  chemical 
source  terms.  Although  it  may  be  inappropriate  to  consider  the 
entire  MILD  combustion  as  a  homogeneous  reactor,  the  elevated 
temperature  of  reactant  mixture  and  distributed  reaction  zones  dis¬ 
cussed  in  Section  3.1  may  generate  local  reaction  zones  having 
B  <  0.  In  the  present  MILD  cases,  not  only  autoigniting  mixtures 
but  also  regions  with  flame  interactions  have  small  C  and  V  yielding 
negative  B  12].  Thus,  regions  with  B  <  0  are  called  reaction  domi¬ 
nated  regions  in  this  study  in  order  to  include  autoignition  and 
interacting  flames.  Clearly,  such  regions  cannot  be  modelled  using 
a  standard  flamelet  approach  involving  scalar  gradient  related 
quantities,  since  the  direct  relation  between  the  scalar  gradient 
and  reaction  rate  no  longer  holds  for  these  and  also  there  is  no  con- 
vective-diffusive-reactive  balance. 

Another  canonical  model  is  a  steady  one-dimensional  premixed 
flame  with  inlet  mixture  velocity  of  SL.  In  this  case,  a  balance 
among  the  convective,  diffusive  and  reactive  terms  gives  B  =  0 
and  the  reactive  contribution  scales  as  prSL/Sth •  In  turbulent  pre¬ 
mixed  flames,  however,  there  are  large  velocity  fluctuations  lead¬ 
ing  to  local  displacement  of  the  flame  and  thus  large  C  and  71  are 
observed  near  flame  fronts  yielding  B  >  0.  This  type  of  reaction 
zone  might  be  similar  to  those  observed  under  Homogeneous 
Charge  Diffusion  Ignition  (HCDI)  conditions  [53]. 

Based  on  these  insights,  the  reaction  and  propagating-flame 
dominated  regions  are  identified  using  B  <  0  and  B  >  0  respec¬ 
tively.  Although,  ignition  front  and  flame  fronts  may  be  identified 
by  computing  displacement  speed  [54]  or  by  using  Chemical 
Explosive  Mode  Analysis  (CEMA)  [55  ,  the  present  analysis 
employs  above  criteria  to  connect  the  identification  more  directly 
with  the  canonical  ignition  and  propagation  phenomena  for  future 
modelling  purposes. 


Figure  12a  shows  typical  contours  of  \C+  -V+\  =  1.5  (dashed 
lines)  and  \7Z+\  =  1.5  (solid  lines)  in  a  y-z  plane  atx+  «  0  (one  grid 
point  inside  the  computational  domain)  for  Case  Bl.  The  contours 
of  | C+  -  V+\  =  1.5  are  wrinkled  and  do  not  seem  to  exactly  follow 
the  contours  of  \7Z+\  =  1.5  although  a  steady  laminar  flame  solu¬ 
tion  is  used  to  generate  the  initial  and  inlet  fields  as  described  in 
Section  2.2.  This  is  due  to  the  role  of  turbulent  mixing  and  molec¬ 
ular  diffusion  that  would  have  occurred  in  the  final  preprocessing 
step  discussed  in  Section  2.2.  However,  these  contours  seem  to 
cluster  at  similar  locations  as  in  Fig.  12a. 

The  comparison  of  these  two  contours  becomes  more  interest¬ 
ing  if  one  considers  a  y-z  plane  at  a  downstream  position.  These 
contours  are  shown  in  Fig.  12b  for  the  y-z  plane  located  at 
x  =  Lx/3.  The  contours  of  \7Z+\  =  1.5  correspond  to  \C+  -  V+\  =  1.5 
contours  in  some  regions  signifying  propagating  flames  locally 
because  of  the  convective-diffusive-reactive  balance.  One  also 
observes  regions  of  predominant  convection-diffusion  and  reac¬ 
tion  alone  without  convection-diffusion.  This  reactive  region  signi¬ 
fies  the  presence  of  autoignition  or  interacting  flames.  The  spatial 
difference  among  these  contours  also  suggests  a  strong  spatial  var¬ 
iation  for  the  unsteady  term  in  Eq.  (11). 

The  local  dominance  of  flame  propagation  or  reaction  is  studied 
by  investigating  the  spatial  variation  of  B+.  The  typical  spatial  var¬ 
iation  is  shown  in  Fig.  13  as  B+  =  1.5  and  -1.5  contours  along  with 
cot .  This  result  is  shown  for  Case  Bl  and  it  is  similar  for  the  other 

Ct 

two  MILD  combustion  cases.  The  reaction  dominated  regions  hav¬ 
ing  B+  <  0  as  noted  earlier  in  this  subsection  are  marked  using  the 
black  contour  lines.  One  of  the  reaction  dominated  regions  having 
B+  <  0  is  marked  using  a  black  box  with  solid  lines,  which  is 
enlarged  for  clarity  (see  Fig.  13,  bottom  right  picture).  As  one  gath¬ 
ers  from  this  enlarged  inset,  B+  <  -1.5  when  co+  is  large.  The 
propagating-flame  dominated  regions,  B+  >  0,  are  marked  using 
the  white  contour  lines  and  one  such  region  is  marked  by  a  white 
box,  which  is  also  enlarged  for  clear  exposition  of  B+  >  1.5  with 
intense  reactions  signifying  the  presence  of  propagating-flames, 
see  Fig.  13,  top  right  picture. 

There  are  many  reaction  zones  having  B+  >  0  for  some  parts  of 
them  and  B+  <  0  for  the  rest.  Four  such  regions  are  marked  using 
black  boxes  with  dashed  lines.  Also,  there  are  regions  with  entan¬ 
gled  reaction-dominated  and  propagating-flame  attributes.  The 
reaction  dominated  events  are  closely  followed  by  flame  propaga¬ 
tion.  The  central  region  marked  using  a  box  with  dash-dotted  line 
in  Fig.  13  is  a  good  example  for  such  a  region  with  both  reaction- 
dominated  and  propagating-flame  phenomena  occurring  close 
by.  This  region  is  also  enlarged  in  the  middle  right  picture  of 
Fig.  13.  The  variation  of  B+  depicted  in  Fig.  13  shows  that  both 
reaction  and  propagating-flame  dominated  regions  coexist  in  MILD 
combustion. 

The  statistical  nature  of  this  coexistence  can  be  investigated 
using  the  PDF  of  B+  and  if  this  coexistence  is  statistically  dominant 
then  this  PDF  is  expected  to  be  uniform.  This  PDF,  typical  for  the 
MILD  combustion  cases  studied  here,  is  shown  in  Fig.  14  for  various 
streamwise  positions  from  the  Case  Bl.  This  PDF  for  the  premixed 
undiluted  case  is  also  shown  in  Fig.  14  for  x+  =  10.2  for  compari¬ 
son.  The  variation  of  the  most  probable  (circles)  and  averaged 
(triangles)  B+  with  x+  is  shown  in  the  inset.  The  behaviour  of 
P(B+)  in  the  premixed  undiluted  case  shows  a  peak  close  to  zero 


dpCT  dpUiCT  _  d 
dt  +  dXi  ~  dXi 

'> - '  s, _ i 
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Fig.  12.  Typical  contours  of  \C+  -V+\  =  1.5  (dashed  line)  and  \TZ+\  =  1.5  (solid)  in  a  y-z  plane.  These  planes  are  extracted  from  an  instantaneous  3D  field  obtained  at 
t  =  1.5td,  and  the  locations  are  (a)  x+  =  0  and  (b)  2.6. 


Fig.  13.  Typical  contours  of  co+  (colour  map)  are  shown  along  with  propagating-flame  dominated  ( B+  =  1 .5,  white  contour  lines)  and  reaction  dominated  ( B+  =  -1.5,  black 
contour  lines)  regions  for  Case  B1  in  the  same  x-y  plane  as  in  Fig.  4d.  A  typical  reaction  dominated  and  propagating-flame  dominated  regions  are  marked  respectively  using  a 
black  box  and  a  white  box  with  solid  lines,  which  are  enlarged  at  bottom  right  and  top  right  respectively.  Several  continuous  reaction  zones  having  both  reaction  and 
propagating-flame  phenomena  are  marked  using  black  boxes  with  dashed  lines.  The  box  with  dash-dotted  lines  in  the  centre  denotes  a  region  where  the  two  phenomena  are 
closely  entangled,  which  is  also  enlarged  at  the  middle  right. 


signifying  the  balance  among  convection,  diffusion  and  reaction  in 
an  average  sense.  The  PDF  is  broader  for  the  MILD  case  and  it 
depends  on  the  streamwise  position,  x+.  The  peak  of  P(B+)  is  close 
to  B+  «  0  for  positions  close  to  the  inlet  boundary.  This  peak  shifts 
to  negative  B+  and  then  moves  towards  B+  =  0  as  the  sampling 
position  moves  downstream.  These  changes  are  gradual  as  shown 
in  the  inset  in  Fig.  14.  The  most  probable  B+,  however,  stays  nega¬ 
tive  as  in  the  inset  but  the  averaged  B+  is  positive  for  various 
streamwise  positions  considered.  Thus,  the  MILD  combustion  is 
dominated  by  strong  unsteady  propagating-flame  behaviour  in 
an  average  sense  although  there  are  many  reaction  dominated 
regions  locally.  The  long  positive  tail  of  P(B+)  resulting  from  con¬ 
vection-diffusion  dominated  regions  yields  B+  >  0.  Similar  behav¬ 


iour  of  this  PDF  is  also  observed  if  the  progress  variable  is  defined 
using  mass  fraction  of  CH4,  H20  or  C02.  The  existence  of  both 
reaction  and  propagating-flame  dominated  regions  along  with  a 
broad  distribution  of  P(B+)  is  also  seen  irrespective  of  the  progress 
variable  definition.  This  signifies  the  coexistence  of  reaction  and 
propagating-flame  dominated  phenomena  along  with  the  impor¬ 
tance  of  strong  local  unsteady  effects  in  the  MILD  combustion. 

It  is  well  known  that  the  gradient  of  cT  plays  an  important  role 
in  turbulent  premixed  combustion  [56-58  .  Thus,  the  scalar  gradi¬ 
ent  magnitude  is  expected  to  be  large  in  propagating-flame  domi¬ 
nated  regions  having  B+  ^  0.  The  contours  of  scalar  gradient 
normalised  using  Sth  are  shown  in  Fig.  15  for  Case  B1  along  with 
the  contours  of  B+  =  1.5  and  -1.5.  The  scalar  gradient  field  is  only 
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Fig.  14.  PDF  of  B+  at  various  x+  positions  in  Case  Bl.  The  PDF  is  also  shown  for  the 
premixed  case  at  x+  =  10.2  (dash-dotted  line).  The  inset  shows  the  variation  of  the 
most  probable  (o)  and  averaged  (A)  values  of  B+  along  x+. 


Fig.  15.  Typical  contours  of  |  VcT|+  (colour  map)  are  shown  along  with  propagating- 
flame  dominated  (B+ =  1.5,  white  contour  lines)  and  reaction  dominated 
( B+  =  -1.5,  black  contour  lines)  regions  for  Case  Bl  in  the  mid  x-y  plane  at 
t  =  1  ,5td.  Regions  with  co+  0 .0  are  masked  using  white  colour  in  the  same  way  as 

in  Fig.  13. 


shown  for  the  regions  with  coCt  ^  1 .0  so  that  a  direct  comparison 
with  Fig.  13  can  be  made.  The  regions  of  high  reaction  rate  are 
not  necessarily  associated  with  regions  of  large  scalar  gradient  in 
MILD  combustion  due  to  interactions  of  reaction  zones  [12]  unlike 
in  premixed  combustion  [57,59  .  This  can  be  seen  by  comparing 
Figs.  13  and  15.  Also,  the  regions  of  large  scalar  gradients  are 
mostly  enclosed  by  white  lines  for  B+  =  1.5  denoting  propagat- 
ing-flame  dominated  regions.  The  substantial  part  of  the  reaction 
dominated  regions  denoted  by  black  lines  in  Fig.  15  have  relatively 
small  scalar  gradient.  These  results  further  reinforce  the  observa¬ 
tion  on  the  distinctive  scalar  gradient  characteristics  of  the  MILD 
combustion  in  a  previous  study  24]. 

In  premixed  combustion,  intense  heat  release  rate  is  associated 
with  large  scalar  dissipation  rate,  which  is  directly  related  to  scalar 
gradient  magnitude  [57-59].  Thus,  regions  with  large  scalar  gradi¬ 
ent  are  conducive  to  the  establishment  of  propagating  flames 
which  propagate  until  they  collide  with  other  locally  propagating 
flames  [12  .  The  propagating  flames  may  not  established  quickly 
in  regions  with  small  scalar  gradient  which  may  be  conducive  to 
establish  reaction-dominated  events  and  thus  these  events  are  pre¬ 
valent  in  these  locations  for  MILD  mixtures.  Similar  observations 
for  reaction  dominated  events  have  been  reported  in  earlier  studies 
of  non-premixed  and  partially-premixed  combustion  [60-62,17]. 

Scalar  gradient  characteristics  in  reaction  and  propagating- 
flame  dominated  regions  can  be  studied  statistically  using  the  PDFs 
of  co+  and  |VcT|+  conditioned  on  B+.  Here,  two  conditions 
B+  ^  -1.5  for  reaction  dominated  and  B+  ^  1.5  for  propagating- 


flame  dominated  regions  are  used.  These  PDFs  are  shown 
in  Fig.  16  for  various  x+  positions  for  Case  Bl.  These  PDFs  are 
constructed  using  all  the  data  sampled  for  this  case. 

The  PDF  of  oj+  conditioned  for  two  conditions  of  B+  is  depicted 

in  Fig.  16a.  It  shows  that  the  range  of  co+  is  broader  for  reaction 

dominated  regions  compared  to  the  propagating-flame  dominated 

regions  and  also  there  is  a  clear  separation  in  the  range  of  oj+ .  The 

PDF  peaks  around  co+  =1.0  for  the  propagating-flame  dominated 

regions  since  oj+  scales  as  prSL/Sth  for  the  premixed  flamelets. 

There  is  a  small  decrease  of  the  value  of  cot  as  one  moves  down¬ 
er 

stream  inside  the  computational  domain.  The  PDFs  are  relatively 
broader  and  the  average  value  of  oj+  is  significantly  larger  for 
the  reaction  dominated  regions.  The  peak  of  the  PDF  decreases  as 
x+  increases,  but  the  range  of  co+  having  non-zero  probability 
increases  suggesting  a  more  uniform  distribution  of  reaction  dom¬ 
inated  regions  in  the  downstream  part  of  the  MILD  combustion 
domain. 

The  conditioned  PDFs  for  |VcT|+  are  shown  in  Fig.  16b.  Although 
the  range  of  |VcT|+  is  similar  for  reaction  and  propagating-flame 
dominated  regions  there  seem  to  be  some  preferred  scalar  gradient 
values  for  these  phenomena.  The  most  probable  scalar  gradient  in 
reaction  dominated  regions  is  generally  smaller  than  that  in  prop¬ 
agating-flame  dominated  regions.  This  offers  further  support  to  the 
physical  picture  gleaned  from  the  DNS  data  discussed  above  for  the 
physics  of  MILD  combustion.  Therefore,  in  the  present  MILD  com¬ 
bustion  conditions,  both  reaction  and  propagating-flame  domi- 


Fig.  16.  PDFs  of  (a)  co+  and  (b)  |Vc7|+  in  propagating-flame  dominated  regions,  B+  ^  1 .5,  and  reaction  dominated  regions,  B+  ^  -1.5,  for  various  x+  positions  given  in  Fig.  14. 
Arrows  indicate  the  general  shift  of  the  PDF  peaks  with  x+.  The  results  are  typical  for  MILD  cases. 
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nated  regions  coexist.  Despite  this,  the  local  scalar  gradient  result¬ 
ing  from  the  mixture  distribution,  turbulence  and  molecular  diffu¬ 
sion  decides  whether  reaction  dominated  or  propagating-flame 
dominated  activities  are  favoured  locally. 


4.  Conclusions 

Three-dimensional  DNS  results  of  turbulent  MILD  and  conven¬ 
tional  premixed  combustion  have  been  studied  to  understand  the 
fundamental  physics  of  reaction  zones  and  their  morphology  in 
MILD  combustion.  The  results  of  the  premixed  case  is  used  for 
comparative  analyses.  The  DNS  configuration  includes  inflow,  out¬ 
flow  and  periodic  boundaries  and  the  inflow  mixture  contains 
pockets  of  reactants  and  products  to  represent  mixtures  for  MILD 
combustion. 

Generally,  MILD  combustion  is  called  flameless  combustion 
implying  that  it  does  not  involve  intense  reaction  zones.  This  study 
using  DNS  data  clearly  demonstrates  that  there  are  regions  with 
strong  chemical  activities  and  these  regions  are  distributed  over 
a  good  portion  of  the  computational  domain  unlike  in  the  pre¬ 
mixed  flame  case  where  the  reaction  zones  are  confined  to  a  small 
portion  of  the  domain.  Also,  interactions  of  reaction  zones  are 
observed  for  the  conditions  of  MILD  combustion  investigated  in 
this  study,  and  the  frequency  and  extent  of  this  interaction  seem 
to  increase  with  dilution.  These  interactions  yield  an  appearance 
of  distributed  turbulent  combustion  under  MILD  conditions  and 
thus  this  combustion  does  seem  flameless. 

The  morphology  of  the  reaction  zones  in  the  MILD  and  pre¬ 
mixed  combustion  is  also  investigated  using  Minkowski  function¬ 
als.  There  are  four  functionals  for  a  three  dimensional  object,  and 
three  characteristic  length  scales  for  the  object  are  deduced  using 
these  four  functionals.  A  shapefinder  map  involving  two  parame¬ 
ters,  planarity  and  filamentarity,  is  obtained  using  the  three  char¬ 
acteristic  length  scales.  These  shapefinders  are  used  for  a 
convenient  and  unambiguous  morphological  classification  of  the 
reaction  zone  structures  observed  in  the  DNS  as  detailed  in 
Section  3.2.  These  structures  are  extracted  using  the  reaction 
rate  value,  corresponding  to  the  maximum  of  Z(£)  = 
(A SNC  cDct\cOct  =  £),  which  is  the  conditionally  averaged  value  of 
the  reaction  rate  weighted  by  the  scalar  dissipation  rate  and  sur¬ 
face  area.  This  avoids  subjectivity  that  may  arise  if  one  uses  a  sim¬ 
ple  threshold  value.  From  this  analysis,  predominantly  sheet-like 
structures  are  observed  for  the  premixed  combustion  case  as  one 
would  expect.  However,  the  MILD  reaction  zones  show  varied 
shapes,  blobs,  pancakes,  small  ribbons,  etc.,  made  of  simply  and 
multiply  connected  objects.  The  statistical  analysis  of  the  shape- 
finders  for  the  MILD  reaction  zones  suggests  that  a  pancake-like 
structure  is  the  most  probable  shape,  which  results  from  events 
such  as  autoignition,  propagating-flames  of  mixtures  containing 
reactants  and  products  with  radicals  and  intermediates,  and  inter¬ 
actions  of  these  flames.  These  phenomena  are  also  entangled  in 
space  and  time  and  separating  them  for  modelling  purposes  may 
not  be  easy.  However,  as  noted  in  Section  3.2.2,  conventional 
turbulent  combustion  models,  such  as  the  EDC  model,  could  be 
improved  by  considering  changes  in  local  reaction  zone  volumes 
due  to  these  phenomena  for  MILD  combustion. 

The  various  fluxes  involved  in  the  progress  variable  transport 
equation  are  also  analysed  for  their  spatial  variations  at  a  given 
instant  and  also  statistically  to  shed  more  light  on  the  above  obser¬ 
vations.  This  analysis  supports  the  above  insight,  suggesting  that 
there  are  regions  dominated  by  reaction  and  propagating-flame 
phenomena  along  with  their  coexistence  for  the  conditions  inves¬ 
tigated  in  this  study.  Although,  these  regions  are  entangled  closely, 
the  local  behaviour  of  scalar  gradient  affects  whether  reaction  or 
propagating-flame  dominated  activities  are  favoured  locally.  These 


activities  produce  reaction  zones  of  varied  shapes  and  sizes  in 
MILD  combustion  and  these  are  analysed  using  Minkowski  func¬ 
tionals  as  noted  above.  These  analyses  show  that  there  are  reaction 
zones  with  holes  resulting  from  contortion  of  reaction  zones  and 
their  interactions  because  of  turbulence.  Local  extinction  cannot 
be  the  cause  for  these  holes  to  form  because  local  extinction  does 
not  occur  for  the  conditions  considered  in  this  study.  These 
interacting  reaction  zones  with  single  or  multiple  holes  offer  new 
challenges  for  turbulent  combustion  modelling.  The  turbulence- 
chemistry  interaction  may  have  some  dependence  on  the  chemical 
kinetic  mechanism  and  thus  the  use  of  a  more  detailed  kinetic 
mechanism  for  the  DNS  is  of  some  interest  for  further  investiga¬ 
tion.  However,  the  broader  perspectives  and  implications  of  the 
observations  made  in  this  study  is  less  likely  to  change. 

The  relative  role  of  autoignition  and  flame-propagation  may 
change  with  a  decrease  in  Tr  and  one  needs  further  direct  simula¬ 
tions  to  assess  this  conclusively.  It  is  also  imperative  that  the  com¬ 
bustion  sub-modelling  for  turbulent  combustion  under  MILD 
conditions  must  be  able  to  cater  for  both  autoignition  and  flame- 
propagation  since  these  processes  are  likely  to  be  interweaved  in 
space  and  time. 
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Appendix  A 

Three  simple  shapes,  a  sphere  with  radius  rs,  a  cylinder  with 
radius  rc,  and  an  elliptic  cylinder  with  re\  and  re2  respectively  as 
minor  and  major  radii  are  tested  to  determine  the  spatial  resolu¬ 
tion  required  to  obtain  accurate  shapefinders.  Figure  A.17  shows 
the  variations  of  P  and  F  for  the  three  test  cases  as  a  function  of 
A Km,  where  A  is  numerical  resolution  (mesh  spacing)  and  Km  is 
the  maximum  curvature  calculated  as  1  /rs  for  the  sphere,  l/rc 
for  the  cylinder,  and  1  /rei  for  the  elliptic  cylinder.  For  large  values 
of  A  Km,  the  fluctuations  in  P  and  F  are  large  due  to  the  lack  of 
spatial  resolution  and  A/cm  ^  0.10  yields  P  and  F  within  ±3.5%  of 
their  actual  (theoretical)  values.  Therefore  in  this  study,  any 
objects  having  A/cm  >  0.10  are  excluded  from  the  calculation. 


Fig.  A.17.  P  and  F  variations  for  test  cases  as  a  function  of  numerical  resolution. 
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